Estimate Hill Numbers
- Bootstrap confidence interval (1000 iterations)
q0 <- iNEXT3D(t(b2), diversity = "TD", q = 0, nboot=1000)
q1 <- iNEXT3D(t(b2), diversity = "TD", q = 1, nboot=1000)
q2 <- iNEXT3D(t(b2), diversity = "TD", q = 2, nboot=1000)
#save(list=c("q0", "q1", "q2"), file="../Data/HillNumbersTD.rda")
load("../Data/HillNumbersTD.rda")
h <- rbind(q0$iNextEst$size_based,
q1$iNextEst$size_based,
q2$iNextEst$size_based)
h <- cbind(h, ta2[match(h$Assemblage, rownames(b2)), 1:2])
h$Month <- factor(h$Month, levels=c("4", "8"))
Fig. X. Size-based diversity accumulation curves based on spe
richness (left), effective number of typical species (middle) and
effective number of dominant species (right). The blue lines indicate
the interpolated (rarefied) and red lines indicate extrapolated parts of
the accumulation curves based on 1000 permutations. Dotted symbols
indicate the observed diversity values.
# Minimum sample coverage for each order q
min(ldply(lapply(splitBy(~Order.q+Assemblage, subset(h, Order.q==0)), FUN=function(x)x[dim(x)[1],]))$SC)
## [1] 0.9957578
min(ldply(lapply(splitBy(~Order.q+Assemblage, subset(h, Order.q==1)), FUN=function(x)x[dim(x)[1],]))$SC)
## [1] 0.9957578
min(ldply(lapply(splitBy(~Order.q+Assemblage, subset(h, Order.q==2)), FUN=function(x)x[dim(x)[1],]))$SC)
## [1] 0.9957578
# SC = 0.99
h1 <- subset(h, Order.q ==0 & SC <= 0.9957578)
p1 <- ggplot()+
geom_ribbon(data=h1,
aes(x=m, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
geom_line(data= subset(h1, Method=="Rarefaction"),
aes(x=m, y=qD, group=Assemblage, colour=Site), size=0.8)+
geom_line(data= subset(h1, Method=="Extrapolation"),
aes(x=m, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
geom_point(data= subset(h1, Method=="Observed"),
aes(x=m, y=qD, group=Assemblage))+
scale_fill_viridis_d()+
scale_color_viridis_d()+
facet_wrap(~Month, scale="free")+
labs(x = "Numbers of Individuals", y = expression(Species~Richness~(""^0*italic(D))),
colour="Site")+
theme_bw()%+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())
h2 <- subset(h, Order.q ==1 & SC <= 0.9957578)
p2 <- ggplot()+
geom_ribbon(data=h2,
aes(x=m, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
geom_line(data= subset(h2, Method=="Rarefaction"),
aes(x=m, y=qD, group=Assemblage, colour=Site), size=0.8)+
geom_line(data= subset(h2, Method=="Extrapolation"),
aes(x=m, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
geom_point(data= subset(h2, Method=="Observed"),
aes(x=m, y=qD, group=Assemblage))+
scale_fill_viridis_d()+
scale_color_viridis_d()+
facet_wrap(~Month, scale="free")+
labs(x = "Numbers of Individuals", y = expression(Species~Richness~(""^1*italic(D))),
colour="Site")+
theme_bw()%+replace% large %+replace% theme(legend.position = "none")
p1/p2

p3 <- ggplot()+
geom_ribbon(data=h1,
aes(x=SC, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
geom_line(data= subset(h1, Method=="Rarefaction"),
aes(x=SC, y=qD, group=Assemblage, colour=Site), size=0.8)+
geom_line(data= subset(h1, Method=="Extrapolation"),
aes(x=SC, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
geom_point(data= subset(h1, Method=="Observed"),
aes(x=SC, y=qD, group=Assemblage))+
scale_fill_viridis_d()+
scale_color_viridis_d()+
facet_wrap(~Month, scale="free")+
labs(x = "Sample Completeness", y = expression(Species~Richness~(""^0*italic(D))),
colour="Site")+
theme_bw()%+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())
p4 <- ggplot()+
geom_ribbon(data=h2,
aes(x=SC, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
geom_line(data= subset(h2, Method=="Rarefaction"),
aes(x=SC, y=qD, group=Assemblage, colour=Site), size=0.8)+
geom_line(data= subset(h2, Method=="Extrapolation"),
aes(x=SC, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
geom_point(data= subset(h2, Method=="Observed"),
aes(x=SC, y=qD, group=Assemblage))+
scale_fill_viridis_d()+
scale_color_viridis_d()+
facet_wrap(~Month, scale="free")+
labs(x = "Sample Completeness", y = expression(Species~Richness~(""^1*italic(D))),
colour="Site")+
theme_bw()%+replace% large %+replace% theme(legend.position = "none")
p3/p4

fr <- strsplit(q0$DataInfo$Assemblage, split="_") %>% ldply
names(fr) <- c("Site", "Month")
fr$Month <- factor(fr$Month, levels=c("4", "8"))
h1_out <- lapply(splitBy(~Assemblage, h1), FUN=function(x)tail(x, n=1)) %>% ldply
h2_out <- lapply(splitBy(~Assemblage, h2), FUN=function(x)tail(x, n=1)) %>% ldply
t0 <- cbind(cbind(fr, S.obs = q0$DataInfo$S.obs), tmp[c(1, 3, 5, 2, 4, 6),])
t0$Site <- factor(t0$Site, levels=c("A", "B", "C"))
t1 <- cbind(h1_out, tmp[c(1, 3, 5, 2, 4, 6),])
t2 <- cbind(h2_out, tmp[c(1, 3, 5, 2, 4, 6),])
h2 <- readRDS("../Data/H2.rds")
h2$Site <- factor(h2$Site, levels=c("A", "B", "C"))
pv <- round(summary(lm(S.obs~Mean, data=t0))$coefficients[2, 4], 3)
r2 <- round(summary(lm(S.obs~Mean, data=t0))$adj.r.squared, 2)
f <- round(summary(lm(S.obs~Mean, data=t0))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)
p7 <- ggplot(data=t0, aes(x=Mean, y=S.obs, shape=Month, colour=Site))+
geom_point(size=2, stroke=1.5)+
stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray", linetype=2)+
annotate(geom = "text", x = 24, y = 12.5, hjust = 0, vjust = 1, label = label, parse = TRUE) +
scale_color_viridis_d()+
scale_shape_manual(values=c(1,19))+
labs(x = expression(Temperature~(degree~C)), y = "Numbers of Species", title="a", shape="Month")+
theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank(), legend.position = "none")
pv <- round(summary(lm(qD~Mean, data=t1))$coefficients[2, 4], 3)
r2 <- round(summary(lm(qD~Mean, data=t1))$adj.r.squared, 2)
f <- round(summary(lm(qD~Mean, data=t1))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)
p8 <- ggplot(data=t1, aes(x=Mean, y=qD, ymin=qD.LCL, ymax=qD.UCL, shape=Month, colour=Site))+
geom_point(size=2, stroke=1.5)+
geom_errorbar()+
stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray", linetype=2)+
annotate(geom = "text", x = 24.5, y = 22, hjust = 0, vjust = 1, label = label, parse = TRUE) +
scale_color_viridis_d()+
scale_shape_manual(values=c(1,19))+
labs(x = expression(Temperature~(degree~C)), y = expression(Species~Richness~(""^0*italic(D))), title="b", shape="Month")+
theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())
pv <- round(summary(lm(qD~Mean, data=t2))$coefficients[2, 4], 3)
r2 <- round(summary(lm(qD~Mean, data=t2))$adj.r.squared, 2)
f <- round(summary(lm(qD~Mean, data=t2))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)
p9 <- ggplot(data=t2, aes(x=Mean, y=qD, ymin=qD.LCL, ymax=qD.UCL, shape=Month, colour=Site))+
geom_point(size=2, stroke=1.5)+
geom_errorbar()+
stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray", linetype=2)+
annotate(geom = "text", x = 24, y = 6, hjust = 0, vjust = 1, label = label, parse = TRUE) +
scale_color_viridis_d()+
scale_shape_manual(values=c(1,19))+
labs(x = expression(Temperature~(degree~C)), y = expression(Shannon~Diversity~(""^1*italic(D))), title="c", shape="Month")+
theme_bw()%+replace% large %+replace% theme(legend.position = "none")
pv <- round(summary(lm(H2~Mean, data=h2))$coefficients[2, 4], 3)
r2 <- round(summary(lm(H2~Mean, data=h2))$adj.r.squared, 2)
f <- round(summary(lm(H2~Mean, data=h2))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)
p10<- ggplot(data=h2, aes(x=Mean, y=H2, ymin=X1, ymax=X2, shape=Month, colour=Site))+
geom_point(size=2, stroke=1.5)+
geom_errorbar()+
stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray")+
annotate(geom = "text", x = 24.5, y = 0.95, hjust = 0, vjust = 1, label = label, parse = TRUE) +
scale_color_viridis_d()+
scale_shape_manual(values=c(1,19))+
labs(x = expression(Temperature~(degree~C)), y = "H2' index", title="d", shape="Month")+
theme_bw()%+replace% large %+replace% theme(legend.position = "none")
(p7+p8)/(p9+p10)

p11 <- ggplot(data=t0, aes(x=Site, y=S.obs, fill=Month))+
geom_bar(stat = "identity", position=position_dodge())+
scale_fill_viridis_d()+
labs(x = "Site", y = "Numbers of Species", title="a", shape="Month")+
theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank(), legend.position = "none")
p12 <- ggplot(data=t1, aes(x=Site, y=qD, ymin=qD.LCL, ymax=qD.UCL, fill=Month))+
geom_bar(stat = "identity", position=position_dodge())+
geom_errorbar(stat = "identity", position=position_dodge(), linetype=2)+
scale_fill_viridis_d()+
labs(x = expression(Temperature~(degree~C)), y = expression(Species~Richness~(""^0*italic(D))), title="b", shape="Month")+
theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())
p13 <- ggplot(data=t2, aes(x=Site, y=qD, ymin=qD.LCL, ymax=qD.UCL, fill=Month))+
geom_bar(stat = "identity", position=position_dodge())+
geom_errorbar(stat = "identity", position=position_dodge(), linetype=2)+
scale_fill_viridis_d()+
labs(x = "Site", y = expression(Shannon~Diversity~(""^1*italic(D))), title="c", shape="Month")+
theme_bw()%+replace% large %+replace% theme(legend.position = "none")
p14<- ggplot(data=h2, aes(x=Site, y=H2, ymin=X1, ymax=X2, fill=Month))+
geom_bar(stat = "identity", position=position_dodge())+
geom_errorbar(stat = "identity", position=position_dodge(), linetype=2)+
scale_fill_viridis_d()+
labs(x = "Site", y = "H2' index", title="d", shape="Month")+
theme_bw()%+replace% large %+replace% theme(legend.position = "none")
(p11+p12)/(p13+p14)

library(knitr)
kable(t0)
| 117 |
A |
4 |
7 |
Hualien Buoy |
24.03 |
121.63 |
4 |
28.0 |
2012 |
25.5 |
22.9 |
2012 |
-1.8 |
3.3 |
0.8 |
| 64 |
B |
4 |
9 |
Chenggong |
23.10 |
121.38 |
4 |
25.7 |
2012 |
24.0 |
22.7 |
2012 |
-2.6 |
0.4 |
-1.3 |
| 32 |
C |
4 |
10 |
Taitung Buoy |
22.72 |
121.14 |
4 |
29.0 |
2012 |
26.6 |
24.9 |
2012 |
-1.0 |
3.1 |
0.7 |
| 121 |
A |
8 |
5 |
Hualien Buoy |
24.03 |
121.63 |
8 |
31.1 |
2012 |
29.2 |
25.3 |
2012 |
-3.2 |
2.6 |
0.7 |
| 68 |
B |
8 |
11 |
Chenggong |
23.10 |
121.38 |
8 |
29.5 |
2012 |
28.1 |
26.9 |
2012 |
-1.6 |
1.0 |
-0.4 |
| 36 |
C |
8 |
8 |
Taitung Buoy |
22.72 |
121.14 |
8 |
30.0 |
2012 |
27.9 |
25.5 |
2012 |
-2.9 |
1.6 |
-0.5 |
kable(t1)
| 117 |
A_4 |
A_4 |
211 |
Rarefaction |
0 |
6.365286 |
5.627048 |
7.103523 |
0.9956038 |
0.9926322 |
0.9985754 |
A |
4 |
Hualien Buoy |
24.03 |
121.63 |
4 |
28.0 |
2012 |
25.5 |
22.9 |
2012 |
-1.8 |
3.3 |
0.8 |
| 64 |
B_4 |
B_4 |
298 |
Extrapolation |
0 |
9.517606 |
7.651172 |
11.384039 |
0.9956365 |
0.9896097 |
1.0000000 |
B |
4 |
Chenggong |
23.10 |
121.38 |
4 |
25.7 |
2012 |
24.0 |
22.7 |
2012 |
-2.6 |
0.4 |
-1.3 |
| 32 |
C_4 |
C_4 |
444 |
Rarefaction |
0 |
7.611364 |
6.480387 |
8.742340 |
0.9954990 |
0.9936961 |
0.9973019 |
C |
4 |
Taitung Buoy |
22.72 |
121.14 |
4 |
29.0 |
2012 |
26.6 |
24.9 |
2012 |
-1.0 |
3.1 |
0.7 |
| 121 |
A_8 |
A_8 |
138 |
Extrapolation |
0 |
5.282364 |
3.289271 |
7.275457 |
0.9956636 |
0.9855805 |
1.0000000 |
A |
8 |
Hualien Buoy |
24.03 |
121.63 |
8 |
31.1 |
2012 |
29.2 |
25.3 |
2012 |
-3.2 |
2.6 |
0.7 |
| 68 |
B_8 |
B_8 |
1428 |
Extrapolation |
0 |
14.932366 |
9.506878 |
20.357854 |
0.9957578 |
0.9924121 |
0.9991035 |
B |
8 |
Chenggong |
23.10 |
121.38 |
8 |
29.5 |
2012 |
28.1 |
26.9 |
2012 |
-1.6 |
1.0 |
-0.4 |
| 36 |
C_8 |
C_8 |
352 |
Rarefaction |
0 |
6.880689 |
5.686621 |
8.074758 |
0.9953146 |
0.9929312 |
0.9976980 |
C |
8 |
Taitung Buoy |
22.72 |
121.14 |
8 |
30.0 |
2012 |
27.9 |
25.5 |
2012 |
-2.9 |
1.6 |
-0.5 |
kable(t2)
| 117 |
A_4 |
A_4 |
211 |
Rarefaction |
1 |
1.883177 |
1.681757 |
2.084597 |
0.9956038 |
0.9925190 |
0.9986886 |
A |
4 |
Hualien Buoy |
24.03 |
121.63 |
4 |
28.0 |
2012 |
25.5 |
22.9 |
2012 |
-1.8 |
3.3 |
0.8 |
| 64 |
B_4 |
B_4 |
298 |
Extrapolation |
1 |
5.129220 |
4.604312 |
5.654129 |
0.9956365 |
0.9896115 |
1.0000000 |
B |
4 |
Chenggong |
23.10 |
121.38 |
4 |
25.7 |
2012 |
24.0 |
22.7 |
2012 |
-2.6 |
0.4 |
-1.3 |
| 32 |
C_4 |
C_4 |
444 |
Rarefaction |
1 |
1.704697 |
1.604944 |
1.804450 |
0.9954990 |
0.9936266 |
0.9973714 |
C |
4 |
Taitung Buoy |
22.72 |
121.14 |
4 |
29.0 |
2012 |
26.6 |
24.9 |
2012 |
-1.0 |
3.1 |
0.7 |
| 121 |
A_8 |
A_8 |
138 |
Extrapolation |
1 |
2.936496 |
2.495292 |
3.377701 |
0.9956636 |
0.9859488 |
1.0000000 |
A |
8 |
Hualien Buoy |
24.03 |
121.63 |
8 |
31.1 |
2012 |
29.2 |
25.3 |
2012 |
-3.2 |
2.6 |
0.7 |
| 68 |
B_8 |
B_8 |
1428 |
Extrapolation |
1 |
3.509970 |
3.259010 |
3.760930 |
0.9957578 |
0.9921777 |
0.9993379 |
B |
8 |
Chenggong |
23.10 |
121.38 |
8 |
29.5 |
2012 |
28.1 |
26.9 |
2012 |
-1.6 |
1.0 |
-0.4 |
| 36 |
C_8 |
C_8 |
352 |
Rarefaction |
1 |
3.880164 |
3.705037 |
4.055290 |
0.9953146 |
0.9928583 |
0.9977709 |
C |
8 |
Taitung Buoy |
22.72 |
121.14 |
8 |
30.0 |
2012 |
27.9 |
25.5 |
2012 |
-2.9 |
1.6 |
-0.5 |